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We describe a model of dynamic Bose-Einstein condensates near a Feshbach resonance that is 
computationally feasible under assumptions of spherical or cylindrical symmetry. Simulations in 
spherical symmetry approximate the experimentally measured time to collapse of an unstably at- 
tractive condensate only when the molecular binding energy in the model is correct, demonstrating 
that the quantum fluctuations and atom-molecule pairing included in the model are the dominant 
mechanisms during collapse. Simulations of condensates with repulsive interactions find some quan- 
titative disagreement, suggesting that pairing and quantum fluctuations are not the only significant 
factors for condensate loss or burst formation. Inclusion of three-body recombination was found 
to be inconsequential in all of our simulations, though we do not consider recent experiments [T] 
conducted at higher densities. 
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I. INTRODUCTION 

Bose-Einstein condensates (BECs) which are sud- 
denly subjected to strong attractive interatomic interac- 
tions, can undergo exotic collapse that resembles super- 
novas . Experiments [51 13] performed in the seemingly 
opposite vein, using strong repulsive interactions, exhibit 
some of the same features of the collapse; namely, an en- 
ergetic burst of atoms, a remnant condensate, and a sig- 
nificant portion of the atoms escaping detection. These 
particular experimental observations of collapsing BECs 
have eluded satisfying quantitative explanation for years, 
with the time to collapse being particularly difficult to 
reproduce in simulations. 

The regime of strong interatomic interactions is 
reached by using a Feshbach resonance, where an exter- 
nal magnetic field allows for tuning the sign and strength 
of the scattering length. We present a theory of BECs 
near a Feshbach resonance, including lowest-order fluc- 
tuations. Since a collapse is a highly local effect, it is 
crucial that we allow for strong inhomogeneities during 
the simulation. This model is defined over time and six 
spatial variables. Symmetry assumptions and the restric- 
tion of our knowledge to only certain off-diagonal corre- 
lations reduces the model to four or five spatial degrees of 
freedom. Assuming spherical symmetry, our simulations 
predict a collapse time of about 2 milliseconds for one of 
the condensates described in Ref. |2j, which agrees well 
with the experiment. Simulations of experiments [3] on 
condensates with repulsive interatomic interactions con- 
sistently overestimate the number of atoms remaining af- 
ter a brief period near the Feshbach resonance. In these 
repulsive simulations, we observe that bound pairs attain 
high velocities (above 8 millimeters per second) immedi- 
ately before dissociating. In all simulations, inclusion of 
an empirically-based model of three-body recombination 



had no significant effect. 

We address the experimentally measured time to col- 
lapse, which has thus far not been accurately reproduced 
for the particular experiments we model. We also explore 
elements of the dynamics that have received little exper- 
imental or theoretical attention, such as the kinetics of 
bound pairs during collapse, and attempt to reproduce 
the results of experiments with a single pulse near the 
resonance. 

Section [ll] describes the experiments and summa- 
rizes past simulations. Section III derives the Hartree- 



Fock-Bogoliubov model and performs the simplifications 
needed to make the resulting equations practical for sim- 
ulation. We present the results of our simulations in Sec- 
tion Hv] and draw conclusions in Section Ivl 



II. OVERVIEW OF EXPERIMENT AND 
THEORY OF COLLAPSE AND RELATED 
DYNAMICS 

A. Experiment 

By exploiting a Feshbach resonance [5], the interac- 
tions between condensed atoms can be tuned from re- 
pulsive to attractive values over only a few microseconds 
[4]. In an often-examined set of experiments [21 El [7] 
conducted at JILA, condensates of about 15,000 ^^Rb 
atoms were formed at a temperature of 3 nanokelvin 
with slightly repulsive interactions. The repulsion was 
balanced by a magnetic trap that was well approximated 
by an axisymmetric harmonic potential, so that an initial 
condensate was stable and neither expanding nor collaps- 
ing. The interactions were then suddenly tuned to be 
attractive. A condensate would appear stable for a short 
time (the collapse time icoiiapsc) after this transition, then 
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suddenly lose atoms at an exponential rate. During col- 
lapse, a burst of energetic atoms was emitted from the 
condensate. Between experiments, the number of atoms 
in the bursts varied by as much as a factor of two, even 
for identical sets of controlled and observed experimental 
parameters. If the atom loss was interrupted by changing 
the strength of the interactions a second time, now to a 
slightly repulsive value, jets of atoms having a lower en- 
ergy than the bursts were emitted, almost entirely in the 
radial direction. A stable, excited, and highly anisotropic 
condensate remained after atom loss ceased. 

A significant number of atoms lost from the condensate 
went undetected; for example, about 8500 atoms out of 
an initial condensate of 15,000 were missing after a col- 
lapse [7. Atoms with energies greater than about 20 ^K, 
atoms in states that were not influenced by the trapping 
potential, and pairs of atoms bound to each other be- 
cause of the Feshbach resonance (henceforth referred to 
as molecules) could not be detected. 

One set of related experiments [3l[7] started from a sta- 
ble, noninteracting or weakly repulsive BEG, which was 
subjected to a rapid magnetic field pulse. The field was 
linearly ramped from its initial value to a value near the 
Feshbach resonance in tens to hundreds of microseconds, 
held at a constant value for one to hundreds of microsec- 
onds, called the hold time, and then quickly and linearly 
ramped back to the initial value. Likewise, the scattering 
length was ramped from zero or a small positive value, to 
a very large positive value, and finally back to its initial 
small value. 

Following a pulse, it was observed that the number of 
atoms remaining in the condensate increased for longer 
ramp times, indicating that the dominant loss mecha- 
nism was not the usual density-dependent loss responsi- 
ble for the rethermalization of a stable condensate. Vary- 
ing the initial density of the BEG did not appreciably 
alter the rate of loss, further suggesting that the loss was 
not density-dependent. As expected, pulses which came 
closer to the resonance resulted in more loss from the 
condensate. When the scattering length was held at a 
large positive value during the hold time, small, damped 
oscillations in atom number were apparent when the hold 
time was varied. 

A burst of atoms similar to that in the collapse ex- 
periments appeared in these experiments with repulsive 
interactions. For small positive values of the scattering 
length, no burst atoms were observed. Varying the num- 
ber of atoms in the surrounding thermal cloud did not 
appreciably affect the bursts, indicating that interactions 
with noncondensed atoms are not responsible. The burst 
atoms remained in the same spin state as the condensed 
atoms, suggesting that spin-flip interactions were not in- 
volved. 

These single-pulse experiments inspired experiments 
m [7] with two magnetic field pulses, separated by a free 
precession time, during which the magnetic field was held 
constant and below the initial value. As with the other 
scenarios, an energetic burst of atoms emanated from 



the condensate. Again, between 8 and 50 percent of the 
atoms escaped detection. 



B. Theory 

BEG collapse has been theoretically studied for sev- 
eral years. Most recently, Altin et al. [T] performed new 
collapse measurements in an optical trap, with a ^^Rb 
condensate that was much denser than those of the JILA 
experiments. In this regime, they found a mean field 
description in combination with three-body losses gave 
a good description of their measured collapse time and 
atom loss curve. In our paper, we focus on the much less 
understood JILA experiments, where three-body losses 
are shown to play an inconsequential role. 

Kagan et al. ^ predicted that collapse occurs on a 
time scale ^collapse ~ where uj is the trap frequency. 
The observations of [5] have shown this prediction to 
be incorrect. Kagan and coworkers also supposed [5] 
that during a collapse, the condensate's density increases 
until density-dependent losses due to three-body recom- 
bination take over, eventually causing expansion of the 
condensate. The cycle then repeats, as the trap pushes 
the remaining condensate back towards the trap cen- 
ter. The GPE simulations of Saito and Ueda and Bao 
PUI [TTl [I3j clearly show such behavior, leading to 
significant atom loss and the prevention of the appear- 
ance of a singularity during collapse. 

These and other [TH [Kl UHl E] simulations qualita- 
tively reproduce the collapse process, the delay before 
atom loss begins, the condensate number decay constant 
Tdocay, bursts, and jets, but have achieved no solid quan- 
titative agreement with observation. Minor differences in 
these authors' results, as well as the lack of quantitative 
agreement with experiment, may be due to their differ- 
ent choices of density-dependent loss rates. These losses 
have been shown [TS' to have a complicated dependence 
on magnetic field, especially near a Feshbach resonance, 
making them difficult to precisely characterize. 

Recognizing the deficiency in atom loss models, Bao 
and coworkers il3i performed a GPE simulation with a 
loss rate chosen so that their simulations correctly re- 
produced the experimental values of ^collapse and con- 
densate remnant number. The atom number decay con- 
stant Tdccay IS reasonably well reproduced, but the simu- 
lated burst energies are much lower than what is exper- 
imentally observed. This failure suggests that a Gross- 
Pitaevskii model with simple density-dependent loss does 
not sufficiently describe the collapse. Savage et al. |19j . 
surveying the literature and performing their own simu- 
lations with several different loss rate coefficients, arrive 
at the same conclusion, noting that theoretical values of 
^collapse are consistently larger than the experimental val- 
ues. The authors mention that this is surprising, since 
the period before collapse begins should be the domain 
where the GPE applies. 

Duine and Stoof '^U\ propose that two condensed 
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atoms can collide, scattering one atom out of the conden- 
sate. They use a Gaussian variational technique to inves- 
tigate this "quantum evaporation" as a possible player in 
the collapse, especially concerning remnant number and 
burst formation. Their simulations show a considerable 
loss from the condensate, but, disagreeing with observa- 
tion, this loss begins immediately after the interatomic 
interactions become attractive. 

Mackie and coworkers [5T] suggest that pairs formed 
by the Feshbach resonance may dissociate into noncon- 
densed atoms during the collapse, and the simulations of 
Milstein et al. [22], which neglect three-body losses but 
include quantum fluctuations and pair formation via the 
Feshbach resonance, show an energetic burst of noncon- 
densed atoms, though using parameters not taken from 
experiments. 

Calzetta et al. |23| downplay the importance of such 
a molecular component for the values of the scattering 
length acoUapsc in the collapse experiments, which are 
far from resonance. Like Yurovsky [21], they attribute 
loss from the condensate to the growth of noncondensed 
modes. Calzetta et al. suggest that a theory account- 
ing for fluctuations would have instabilities growing out 
of those fluctuations, which may account for earlier col- 
lapse times, or that a loss of coherence between atoms is 
largely responsible for atom loss [IS]. 

Wiister and coworkers (25] use the same theory of fluc- 
tuations as Milstein et al. but also regard the molecular 
component as unimportant to collapse. Their simulations 
still find a tcoiiapse exceeding the observed value. Using 
an alternate, truncated Wigner formulation and includ- 
ing initial and dynamical noise, a background thermal 
component, and cylindrical geometry with experimental 
parameters, Wiister et al. [22] still overestimate the ex- 
perimentally measured tcoiiapsc by about 40 percent. 

Haldar et al. 28 summarize a correlated potential har- 
monics expansion method that accounts for two-body 
correlations and models interatomic interactions with the 
van der Waals potential. They show that anharmonicity 
and a finite potential barrier at the ends of the optical 
trap have a nontrivial effect on the stability of attractive 
condensates. The same method is used to demonstrate 
a variation of energies with effective scattering length 
where mean field theory predicts none |29| . and predicts 
the critical number of condensed atoms at which a con- 
densate collapses much more accurately than does mean 
field theory ^30j, highlighting the importance beyond- 
mean-field effects in collapsing condensates. 

All these models have at least some qualitative agree- 
ment with observation, and some provide insight into 
other aspects of the collapse experiments [12 [T31 ITS] . 
Saito and Ueda [10] suggest the bursts are atoms origi- 
nally near the center of the collapse that acquire kinetic 
energy when three-body losses suddenly remove a large 
number of atoms from the center of the collapse. In these 
simulations and others [HI US HH HH [HI [l^ , the burst 
atoms are distinguished from the condensate purely by 
their location. In the simulations of Milstein et al. and 



Wiister et al. the burst is assumed to be a distinct non- 
condensed field, which can occupy the same space as the 
condensate. 

The magnetic field pulse experiments have stimu- 
lated fewer simulations than the collapse experiments. 
Duinc and Stoof [3T] use coupled mean fields allow- 
ing for molecule formation, quantum evaporation, and 
three-body losses in modeling the one-pulse experiments. 
These simulations had only general qualitative agree- 
ment with the experiments, but with the observation that 
the inclusion of three-body losses supressed oscillations 
in numbers of atoms and molecules, despite the belief 
that these density-dependent losses should be unimpor- 
tant under the experimental circumstances [7]. Mackie 
and coworkers [2J use a coupled mean field model that 
allows for dissociation of molecules into noncondensed 
atom pairs, but find only about five percent loss to the 
noncondensed component, with very few molecules be- 
ing retained. The authors observe a larger loss in simu- 
lations of the two-pulse experiments, but the oscillation 
envelopes have a behavior markedly different from the 
slight damping observed in the experiments. Kokkclmans 
and Holland [32] use the same model as the Milstein et al. 
collapse simulation [55] , but use a Gaussian average over 
a homogeneous gas to simulate the behavior of a trapped 
gas. These simulations agreed fairly well with the two- 
pulse experiments, showing that the majority of atoms 
lost from the condensate go into noncondensed modes, 
and the missing atoms are identified as molecules. Kohler 
and coworkers |33j model the two-pulse experiments with 
a theory that includes molecule formation and quantum 
fluctuations, and find good qualitative agreement with 
the experiments. They also find that the presence of the 
trap moves the means of the condensate and burst num- 
bers' oscillations closer together in a way not captured 
by a Gaussian average of a homogeneous gas. They at- 
tribute this difference to the presence of a length scale 
not found in the homogeneous gas simulations. 

Many of the questions posed by the experiments have 
eluded satisfying explanation. In the case of collapse, 
the experimentally measured tcoUapso has been particu- 
larly difficult to simulate. Unsettled points of contention 
include the mechanisms by which the jets and bursts op- 
erate, and the importance of three-body losses to the col- 
lapse. The counter-intuitive behavior of the condensate 
after the collapse has thus far received relatively little 
attention |34| . as have the simulations involving a sin- 
gle pulse of repulsive interactions. There has also been 
little exploration of the various models' implications for 
measurable quantities and phenomena that have not yet 
been vigorously pursued in experiments. 



III. HARTREE-FOCK-BOGOLIUBOV MODEL 

To treat both collapse experiments and single and dou- 
ble pulse purely repulsive experiments we will work with 
a Hartree-Fock-Bogoliubov (HFB) model which explicity 
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takes into account the two main channels involved with 
the Feshbach resonance. Previous versions of this model 
exist [221 |3S], including operator equations and factor- 
ized expectation values of the HFB equations f551 
but we present adaptations to spherical and cylindrical 
geometry in a complete and rigorous derivation, as not 
found elsewhere to our knowledge. 

We begin with a definition of the atomic field operator 

as 



V'q(x) = 0a (x) +X(x) 



(1) 



where (/>a(x) = (V'a(x)), and all operators are taken to be 
in the Heisenberg picture, though for brevity we do not 
explicitly write the time dependence. We call a pair of 
atoms in a quasibound state due to the Feshbach reso- 
nance a molecule. We loosely assign the term "molecule" 
to the closed channel for the purpose of this discussion. 
The actual molecule is a superposition of open and closed 
channels. The field operator for molecules is represented 
by ■(/'m(x), with expectation value (/'m(x) = (^/'^(x)). 
To keep the problem tractable, we assume no fluctua- 
tions around the molecular condensate. The normal and 
anomalous fluctuations are deflned by 



G'Ar(x,x') 



and 



Ga(x,x') 



(X^(x')x(x)) 

(Vil(x')V'a(x)) 

(X(x')x(x)) 



(V31(x'))(^.(x)) 

0:(x')0a(x) (2) 



(^,(x'))(V3a(x)) 
0,(x')0a(x) (3) 



respectively. The diagonal elements (x = x') of the nor- 
mal fluctuations give a number density of noncondensed 
atoms, while the diagonal elements of the anomalous fluc- 
tuations give the variance in the mean 0a(x). Off diag- 
onal elements of the normal and anomalous fluctuations 
are equal-time correlation functions. 

We obtain equations of motion for the atomic and 
molecular mean fields and the normal and anomalous 
fiuctuations by finding Heisenberg equations of motion 
for the atomic and molecular field operators and for the 
products x^(x')x(x) and x(x')x(x). Taking the expec- 
tation value on both sides of the resulting equations re- 
sults in averages of products of atomic and molecular 
field operators; by assuming that the atomic and molec- 
ular field operators act in orthogonal subspaces of the full 
Hilbert space, these expectation values factorize as, for 
example, 

(^t(x)v;„,(x'))-0:(x)0,„(x'). (4) 

Products of three or four atomic field operators, which 
do appear in the equations of motion, require special care 



to be expressed in terms only of the atomic mean field 
and normal and anomalous fiuctuations. 

We may apply a manifestation of Wick's theorem [37] 
if the state of the system is an eigenstate of every Bogoli- 
ubov atomic quasiparticle annihilation operator (a linear 
superposition of the atomic momentum-space creation 
and annihilation operators). This quasiparticle coherent 
state, which is not generally a coherent state of the field 
operator ipa{^), is a squeezed state [35], and thus satu- 
rates the number-phase Heisenberg uncertainty relation. 
Then 



(x) V-, (x) (x) ) = I (x) r 0„ (x) + (x) (x, x) 



+ 20„(x) GAr(x,x) 



(5) 



for example, is exact. 
We use the Hamiltonian 



H = 



d^^"V^l(x") 
d^:^'>^„(x") 



--V^ + nx") 



V3a(x") 



-^v2 + 2y(x") + i^ 

4m 



v;^„(x") 



u 



d3a;'>t(x")^t(x")v;,(x")V'a(x") 



^1„(x")^a(x")V'a(x")+h.C. , (6) 



which is often referred to as a two-channel Hamiltonian, 
where m is the mass of an atom, V{tc") is the exter- 
nal potential felt by a single atom, v is the detuning of 
the Feshbach resonance, U relates to the strength of the 
non-resonant atom-atom interaction, and g relates to the 
strength of the atom-molecule coupling which gives rise 
to the Feshbach resonance. These interaction parameters 
are based on the assumption of contact interactions be- 
tween the particles. Such interactions, however, give rise 
to an ultraviolet divergence in momentum space, which 
must be treated properly by renormalization. This is 
done by the introduction of a momentum cutoff K, while 
making sure that at the same time the correct underlying 
two-body resonance physics is maintained. Therefore it is 
necessary to consider the contact interactions as the zero 
range limits of the actual nonlocal interactomic poten- 
tials. The properties of the contact potentials can then be 
chosen such that the two-body physics around a Feshbach 
resonance is correctly described [32l [36] . This renormal- 
ization procedure amounts to a if -dependent relationship 
between the interaction parameters in the Hamiltonian 
and the physical interaction parameters given by 

U = FC/o, 
9 = Tgo, 



where 



1 



1 - aUo 



and a ; 



mK 
2^' 



(7) 



(8) 
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Parameter 


Value 


m 


84.91 1794mproton 


abg 


-450.0ao 


AS 


10.95 G 




-2.2259^is 




154.9 G 



TABLE I: Values of the fixed renormalization parameters that 
we used in our code, where m^rotou is the mass of a proton, 
ao is the Bohr radius, and /is is the Bohr magneton. 



The parameters with a subscript are the unrenormal- 
ized physical interaction parameters and are defined by 



m 



go = VgcUo AB An 



(9) 



where abg is the background scattering length; AB is 
the width of the Feshbach resonance, defined as the dis- 
tance from the resonance position to the point where the 
effective scattering length is zero; A^^^^ is the differ- 
ence in magnetic moments between an uncoupled bound 
and unbound pair of atoms; B is the external magnetic 
field; -Bics is the position of the resonance, defined as 
the value of the magnetic field for which the effective 
scattering length diverges; and the correction factor gc 
may be set to 1.816 to match the binding energy of the 
contact potential model as closely as possible to the field- 
dependent binding energy of the weakly-bound rubidium 
Feshbach state jS^, or to 2 to match the mean field en- 
ergy. The values of those parameters which are fixed are 
summarized in Table |l| The cutoff K is set to 6 x 10* 
inverse meters, the largest wavenumber in our simula- 
tions. One may combine the unrenormalized parameters 
^ and floff = Obg [1 - AB/ [B - B^^^)] to write 



9l 



(10) 



Consistent with the two-body approximation, the 
Hamiltonian (|6| neglects all interactions between atoms 
and molecules that do not involve molecule formation or 
dissociation, and assumes molecules do not interact with 
each other. 

When we transform to the center of mass R = 
(x + x') /2 and relative r = x — x' coordinates, we can 
write 

4(R) = 0,(x) = </.,(x'), 

0™(R) = 0m (x) = 0„i(x') , 

Gn{^,v) = Gn (x,x'), 
GA(R,r) = G^ (x,x'), 

y(R)^nx) = nx')- (11) 



These statements will be valid if, in the case of spheri- 
cal geometry, the external potential and the initial con- 
ditions on all single-particle fields (diagonal elements of 
Gn and Ga included) are rotationally invariant, and we 
only consider x and x' such that 

|x| = |x'| . (12) 

For cylindrical geometry, we assume the initial conditions 
on single-particle fields and the external potential are in- 
variant with respect to rotation about a vertical axis (let 
it be the z-axis) and invariant with respect to reflections 
over a plane normal to that axis (the x-y plane); we must 
also restrict x and x' such that 

l^pl = I^pI : 

|x.i-|x;i, (13) 

where Xp is the component of x lying in the x-y plane 
and yiz is the component along the z-axis. Note that 



the restrictions ( 12 ) and ( 13 1 impose no additional ap- 



proximations beyond those that have already been made. 
They merely provide convenient simplifications in the 
equations of motion. 



The notation in Eq. (Ill is general enough to handle 
the spherical and cylindrical cases, though only the mag- 
nitude of R is important to single-particle fields in the 
former, and only the magnitudes of Rz and Rp are impor- 
tant to single-particle fields in the latter. For consistency, 
we make the same assumptions about the dependence of 
off-diagonal correlations on R for each symmetry. Fol- 
lowing |22j . we then Fourier transform over the relative 
coordinate: 



G^(R,r)^Gjv(R,k), 
GA(R,r)^GA(R,k), 



(14) 



which are valid because Gn and Ga are each symmetric 
with respect to r. This transform removes a Dirac delta 
function appearing in the partial differential equation for 
the anomalous fluctuations. 

Next, we choose coordinate systems appropriate to the 
geometry. In the spherically symmetric case, we choose 
spherical coordinates, where k is aligned with the 2-axis, 
as in Figure [l] The assumption of spherical symmetry in 
R permits this choice, since only the relative orientation 
of R and k will be important. The angle between R and 
k is then 0, and the azimuthal angle of R in this coordi- 
nate system is (j). For cylindrical symmetry, we are only 
permitted to rotate the coordinate axes about the z-axis 
without changing the values of each dependent variable, 
so we align the z-axis with kp, the component of k ly- 
ing in the x-y plane, as in Figure [2] The flve spatial 
variables are then k^, kp, Rz, Rp, and 0, the axial and 
radial components of the relative wavenumber and center 
of mass coordinate, repectively, and the azimuthal angle 
of R. The Laplacians and gradients involved are then ex- 
pressed in spherical or cylindrical coordinates, depending 
on the geometry; in spherical symmetry, the radial parts 
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FIG. 1: Coordinate axes for spherical symmetry. Since all 
fields are independent of the orientation of R, we are free to 
rotate the axes (keeping the origin fixed), and so we align 
the z-axis with k. Then the dependence of the correlation 
functions' values on the relative orientations of k and R (bold 
arrows) can be expressed in spherical coordinates. 




y 



FIG. 2: Coordinate axes for cylindrical symmetry. Since all 
fields are independent of the azimuthal angle of R, we are free 
to rotate the axes so long as the 2-axis and origin remain fixed. 
We align the i-axis with the component of k lying in the x- 
y plane. Then the dependence of the correlation functions' 
values on the relative orientations of k and R (bold arrows) 
can be expressed in cylindrical coordinates. 



of the Laplacians may be simplified, as in the usual treat- 
ment of the hydrogen atom, with the substitutions 

gNiR:k,e,^) = RGNiR,k,9,^), (15) 

and likewise for the molecular field and anomalous fluc- 
tuations. In cylindrical symmetry, we use 

Va{Rz,Rp) = y/Rp^t'aiRz, Rp) , 

GN{Rz,Rp,kz.,kp,4)) = y/RpGN{Rz,Rp,k^,kp,4)) , 

(16) 

and likewise for the molecular field and anomalous fluc- 
tuations. 

Note that, in the spherical case, we have not assumed 
that the azimuthal angle is unimportant, making our 
model more general than those used earlier [521 
Rather than the Legendre polynomial expansions used by 
Milstein et al. [22] and Wiister et al. [3S] , an appropriate 
basis for expanding the normal and anomalous fluctua- 
tions is then the spherical harmonics Y^'^[6,(j3). We write 



oo I 

GN/AiR, k,e,cj,) = J2Yl ^N/AiR^ ^) ^;'(^' ' (17) 

1=0 q=-l 



where the N/A subscript means that the equation applies 
to both Qisi and Qa- Only one angle is present in the 
cylindrical case. Therefore, we use trigonometric func- 
tions in 4> to form a complete angular basis, which is com- 
mon for spectral solutions to PDEs [35]. In this respect, 
sine series are slightly more stringent in their criteria for 
uniform convergence than cosine series. Accordingly we 
expand the normal and anomalous fluctuations as 



Gn/a{Rz,Rp, kz,kp, (j)) = 

OO 

Y,GN/AiRz^Rp^k,,kp) cos {n(f>) (18) 

n=0 



where a superscript n indexes a generally complex scalar 
expansion coefficient and does not denote a power. 

The HFB model for spherical symmetry in the center 
of mass coordinate consists of the dcnumerably infinite 
set of equations 
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d 
d 



d 

"at' 



+ Va{R) + U 



E? 

iPm{R) + 



m 



{l-q + l){l + q+l) 
{21 + 1) {21 + 3) 



R 

vl{R) 

R 
l + l 



MR) ■ 



^rSiiR) I g^m{R) 



R 



d_ 

dR'^ R 



{I + q)il- q) 



'+iY (2/ + l)(2Z- 1) \dR R 



+ Ga{R) 



Q'^^'\R,k) 



gN'''{R,k) 



MR) , OiiR) 
R^ R 



Q*^'-\R,k) 



'<\R) , Q*a'{R) 



ih^/i\R,k) 



+ 



R 

_ 

4rn 

<fm{R) 

' R 



+ u 



J?2 

I {l + l) 
i?2 



+ 



R 



Q'fiR^k), 
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g'/{R,k) 



g'^\R, k) + {-ly g*J'-''{R, k) + v^^o.z 5o,qR 



(21) 



(22) 



where the diagonal elements ^/^/^ of the normal and 
anomalous fluctuations, respectively, are computed from 

1 

^n/a{R) = dA: k'g°//AR^ k) . (23) 



The discrete step function u„ is 1 if a > 6 and otherwise. 
In the cylindrical case, the denumerably infinite set of 
equations is 
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Qn/a{Rp^Rz) 



(27r)^ 

/•OO 

/ kp dkp Q j^,j^[Rp, Rz, kp, kz) 
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Note that the normal fluctuations will be real in either 
symmetry if they are initially real. 

Where desired, we can incorporate a three-body loss 
coefficient into our spherical model for the atoms inside 
the atomic condensate, by fitting the analytical model 



used by Braaten and Hammer [40] and D'Incao et al. [41] 
to the empirical data of Roberts et al. [18] . This model 
assumes that universal Efimov physics, in combination 
with the weakly-bound Feshbach state and deeply bound 
dimer states, give rise to the three-body recombination. 
The result is the addition of the term 



MR) 



to the right hand side of Eq. (19), where 



(28) 



K3 (floff) 



h 4 Jer.l • e-2'' [cos^ (solnacff/flpos) + sinh^T?] + 16.8 (1 - e 
2m I 4590 (sinh 2ri) / [sin^ (sq In flcff/ancg) + sinh^ 77] 



if Ocff > 
if Ocff < 



(29) 



with So — 1.00624, ry = 0.0165971, Opos = 236.743ao, and a„og = — apos/0.96. Figure [s] compares this loss rate with 
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FIG. 3: A comparison of the measured [TSl (circles) and ana- 
lytical (solid line) three-body loss rate of Eq. ( 29 1 as a function 
of magnetic field. Eq. ( 29 1 has been multiplied by a factor of 



6 here, because the measurements were done in a thermal gas. 
The dashed vertical line marks B = 154.9 Gauss, where the 
effective scattering length diverges. 



measured values [18]. The leading scaling goes as a|Jg, 
which has been shown to be effective elsewhere [H]. In 
practice, we impose an empirically based floor of 10~^* 
cm^/s on K3{acs). 



IV. RESULTS OF HFB SIMULATIONS 

Our simulations are performed in spherical symmetry 
using the method of lines [43] . We do not perform more 
realistic axisymmetric simulations, since they would re- 
quire a larger amount of time. The infinite sums over I 
in Eqs. 



21) and (22 1 are truncated at ^ = 1, which is 



the smallest value that provides satisfactorily converged 
results. Truncating at Z = 3 docs not significantly affect 
the results (for example, atom number during a collapse 
simulation changes by less than 0.002 percent and by 
less than 0.003 percent during the longest, most repul- 
sive pulse simulation), though computation time is signif- 
icantly larger. The expansion coefficients for any nonzero 
q are coupled only to those for — q, and only g = is 
necessary for calculation of the diagonal elements of the 
fluctuations. Therefore, we only propagate q — in our 
simulation, and forego whatever information about fluc- 
tuations is present in the q coefficients. A fifth order 
Runge-Kutta method provides time propagation, and we 
handle the spatial derivatives by spectral collocation in 
a sinusoidal basis. We neglect three-body effects unless 
noted otherwise: as we will show, three-body effects do 
not significantly alter our results. 

We create initial states (non-interacting in attractive 
simulations and with = 7ao in repulsive simula- 



tions) using imaginary time propagation with the Gross- 
Pit aevskii equation in the trap used in Ref. using 
w = 27r X 12.77 radians per second, which is the geomet- 
ric mean of the experimental trapping frequencies. Imag- 
inary time propagation is necessary because the simula- 
tion's implementation imposes infinite potential barriers 
at the end of the grid, so that a closed-form solution for 
the Gross-Pitaevskii equation is not available. We as- 
sume the molecular field's initial state is 



2^' 



(30) 



which is the exact initial state for a uniform gas with- 
out fluctuations. Noting the importance of the molecular 
binding energy [32) . we always choose gc = 1.816, except 
where noted. 



A. Attractive Interactions 

The collapse simulation begins with the effective scat- 
tering length set to — 12ao, where it remains for the du- 
ration of the simulation. All simulation parameters are 
taken as the physical parameters of the experiment of 
Ref. [2]. 

Figure [4| plots the number of atoms in the condensate 
as a function of time, and Figure [5] shows the conden- 
sate density. The erratic oscillations around the gen- 
eral downward trend visible in Figure [4] are more clearly 
shown in Figure [6] We interpret 2 milliseconds as our 
predicted collapse time, because the condensate num- 
ber begins to drop substantially and the density grows 
increasingly concentrated near the origin at about this 
time. This value can be directly compared to the nearest 
data point in Donley et a/.'s Figure 2 [5], which has a col- 
lapse time of about 2 ± 1 milliseconds. Beyond 2.04 mil- 
liseconds, negative noncondensed densities consistently 
start to appear in our simulation, and the error in total 
number diverges. Such instability may occur because the 
spatial grid is uniform and non-adaptive, while most of 
the dynamics occur on the innermost grid points once 
the collapse has begun. Thus, an adaptive spatial grid 
would be required to be sure of a ^coiiapso value more 
precise than 2 milliseconds. Including three-body effects 
changes the results by less that 0.2 percent, as Figure [7] 
shows. The K3 factor in the collapses is at the floor value 
of 10"^^ cm^/s, which is twice the maximum value that 
Altin et al. [T] find allowable. However, since the loss 
term scales as the square of the density, and Altin et al. 's 
condensates are several times denser than those of the 
JILA experiments, it seems likely that would indeed 
have to be much larger to have a noteworthy effect in our 
simulations. 

These results are compared to mean field theoretic re- 
sults in Figure [Sj which plots atomic condensate density 
at the origin as calculated from simulations using our 
model and the GPE. A comparison of numbers of atoms 
would be pointless, since the GPE only allows particles 
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to be in the atomic condensate. Characterizing collapse 
by density at the origin, the GPE predicts a collapse time 
more than 5 milliseconds later than our model predicts 
and the experiments measure. Figure [8] also shows the 
results obtained with our model when B — Sics + AB. 
First, we use the nominal value of the coupling correction 
factor 5c = 2. According to Equation (10 1, the scatter- 
ing length Qcs = for this field. Then, both the GPE 
simulation and our model simulation result in a station- 
ary solution, as would be expected for the noninteracting 
case. When we set the scattering length to — 12ao, our 
model gives results nearly identical to the GPE on atom 
number and condensate density at the origin (see Fig- 
ure |8]), and a collapse occurs in both simulations. These 
results are consistent with the findings in Ref. |26J . How- 
ever, when we tune the coupling correction factor to the 
value gc = 1.816 to match the correct molecular binding 
energy, the simulation results of our model differ signifi- 
cantly from the nominal GPE results. Simulation results 
are again shown in Figure [8]for B = Brcs + A-B, but now 
we observe a collapse. This can be understood again from 
Equation ( 10 1, since at this field for = 1.816 we obtain 
the value Oes = —AIuq. 

As we discussed before, it is not possible in our model 
to have the binding energy and the mean field energy 
both correctly described, which leads to a trade-off in 
what is the dominant effect causing the collapse. For our 
treatment of the collapse problem, we argue that the col- 
lapse near a Feshbach resonance is most sensitive to the 
binding energy; however, since we have to set accordingly 
the value gc — 1.816, this unfortunately leads to stabil- 
ity problems in the zero-crossing region of the scattering 
length. 

The molecular field is very sparsely populated, never 
totalling more that 10 percent of a molecule, even though 
pair formation is the primary mechanism in the model 
[36j of interatomic interactions that we used. However, 
this does not mean the molecular field could simply be 
neglected: in a two-channel model of a Feshbach reso- 
nance the molecular field often has very small occupation 
but plays a key role in both static and dynamic descrip- 
tions of the physical molecule, a superposition between 
the atomic (open) and molecular (closed) channels. The 
oscillations between atoms and molecules, especially vis- 
ible in Figure [6j have a period of 1.1 microseconds, which 
is approximately equal to tjj = m?U'^ /fv' , a natural time 
scale arising from dimensional analysis of the parameters 
of the model. Other analytically deduced time scales are 
tg = h7 /m^g'^ « 0.57 nanoseconds, and = h/i' « 22 
nanoseconds. 

While most of our collapse results are for Ocfr — — 12ao, 
Figure [9] shows atomic condensate density at the origin 
for other values of Ucg. We find that for a^tf — — 26.5ao 
and floff = — 54.5ao, respectively, the simulations become 
unstable after 1.88 milliseconds and 1.6 milliseconds, re- 
spectively. These times are within the error bounds for 
the corresponding experimental collapse times reported 
in Figure 2 of Ref. [21 , and correctly capture the general 
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FIG. 4: Number of condensed atoms in the collapse simula- 
tion. The dashed vertical line marks t — 2.04 ms, the last 
sample point for which we simulation is stable. 
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FIG. 5: (Color online) Condensed atom density in the collapse 
simulation. In the simulation, R extends to 10 /im. 



trend of collapse time decreasing nontrivially as effective 
scattering length becomes more negative. 



B. Repulsive Interactions 

It is interesting to compare our nonhomogeneous 
model with the double pulse experiment of Donley et 
al. 4 , which has been described succesfully earlier with a 
similar but homogeneous model [35] by applying the local 
density approximation. In this experiment, a condensate 
is subjected to two magnetic field pulses with interatomic 
interactions that become very repulsive. Although we do 
not expect as large inhomogeneity effects as in the at- 



11 




t (ms) t (ms) 



FIG. 6: Number of condensed atoms during a brief period in 
the middle of the collapse simulation. The plotted points are 
0.5 nanoseconds apart, so the visible oscillations of period 1.1 
microseconds are not aliased. 
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FIG. 7: Difference over the average of the number of con- 
densed atoms between collapse simulations with and without 
three-body effects. The dashed vertical line marks t = 2.04 
ms, the last sample point for which the simulation is stable. 
At no time does the average difference exceed 0.2 percent, 
though it is consistently increasing. 



tractive situation where a collapse occurs, it is useful to 
compare the results between the models as it indicates 
the validity of the local density approximation. The re- 
sults are indeed different. With an initial state having 
16, 600 atoms, all in the condensate, and a free evolution 
time of 10 microseconds (the left-most set of data points 
in Figure 5 of Ref. [31]), we find about 12, 000 condensed 
atoms and slightly over 4, 500 noncondensed atoms at 



FIG. 8: Density of atoms in the atomic condensate as pre- 
dicted by the Gross-Pitaevskii equation (solid line) during a 
simulated collapse and by various configurations of our model. 
For readability, the abscissa here is truncated before the GPE 
simulation ends; its density at the origin climbs until insta- 
bility terminated the simulation at 7.35 milliseconds. Our 
model with Qc = 1.816 (long dashed line) collapses much ear- 
lier than the GPE for a given magnetic field, but also collapses 
for B = Bros + A_B (short dashed line). For Qc = 2, our model 
predicts a rise in density very similar to the GPE (a plot 
would overlap the GPE result shown here), and no effective 
change in density when B — Bros + AS (dotted line). The 
"C" and "S" labels in the legend are reminders as to which 
series are for nominally collapsing and stationary values of 
flcff, respectively. 
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FIG. 9: Atomic condensate density at the origin for simula- 
tions with effective scattering lengths of (solid line), — 12ao 
(long dashes), — 26.5ao (short dashes), and — 54.5ao (dotted 
line). These simulations become unstable after 2.31, 2.04, 
1.88, and 1.6 milliseconds, respectively. All simulations use 
gc = 1.816. 
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FIG. 10: Number of condensed atoms in the two-pulse simu- 
lation. The oscillations in the middle of the simulation occur 
during the free precession time after the first pulse and before 
the second, when the magnetic field is held constant. 



the end of the simulation. Kokkelmans and Holland |32j 
find about 9,500 and 6,300, respectively, while the ex- 
periment [4l ends with just over 10,000 condensed atoms 
and slightly under 5, 000 noncondensed atoms. Figure [T0| 
shows the number of condensed atoms over the course of 
our simulation. It should be noted that a homogeneous 
version of our code exactly reproduced the results of the 
code used in Ref. [52] . 

In comparison with a second type of experiment, we 
conducted simulations using the same parameters as 
those used to create the left-most three sets of data points 
in Figure 4 of Ref. [3], in which a condensate is subjected 
to a single brief magnetic field pulse of strength Bp^ise- 
Both the pulse strength and the time iramp to reach that 
value of magnetic field are varied between experiments. 



Figure 11 shows the time dependence of the magnetic 
field used in the simulations. The results are plotted 



alongside Claussen et aVs ^Mj in Figure 12 which in- 
cludes experimental error bars. Error bars on the sim- 
ulated data points, which would reflect the differences 
between different grid resolutions, are too small to see in 
the figure. Our simulations consistently end with higher 
numbers of condensed atoms than the experiments did, 
and the condensed atom number decreased monotoni- 
cally with longer ramp times. However, the experimen- 
tal results are not easy to interprete. For instance, there 
is a local minimum in the 158 Gauss series which is not 
observed in the simulations. Moreover, the experimen- 
tal data have a larger number of atoms remaining for 
the 156.7 Gauss pulse than the 157.2 Gauss pulse for the 
shortest ramp time; our simulations do not reproduce this 
effect, either. Including three-body effects does not sig- 
nificantly change our results, as Figure [T3| shows. Again, 
the term in these simulations is larger than what Al- 
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FIG. 11: Time dependence of the magnetic field for simula- 
tions with a single pulse that induces repulsive interatomic 
interactions. The lowest value of the magnetic field is Bpuisc, 
and the time to descend from the initial value Bini ~ 165.68 



Gauss, where Ooff = 7ao, to -Bpuise is called 



Our simula- 



tions had the magnetic field near resonance for 1 microsecond, 
though experiments were conducted j3] with many other val- 
ues. 
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FIG. 12: Numbers of condensed atoms remaining at the ends 
of simulations (glyphs connected by dashed lines) and exper- 
iments O [44] (glyphs connected by solid lines) as a function 
of tramp. Curvcs are provided as a guide to the eye; actual 
data is represented by glyphs. 



tin et al. predict, but note that these simulations are for 
relatively short times, and for condensates that are over 
an order of magnitude sparser than those of Altin et al. 
and expanding further. 

Our model's results for atomic condensate density at 
the origin are compared to those of the GPE in Figure [T4| 
for the -Bpuiso = 156 Gauss case. Not only does the GPE 
predict an invariant number of condensed atoms, but it 
also predicts a nearly invariant density. This comparison 
is representative of all of our simulations with repulsive 
interactions. 

Figure [15] shows the molecular velocities during a small 
portion of the i?puise = 156 Gauss simulation, calculated 
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FIG. 13: Difference over the average of the number of con- 
densed atoms between simulations with and without three- 
body effects. The scenario was a single pulse to -B = 156 
Gauss. At no time does the average difference exceed 2.7 per- 
cent, though it increases most dramatically when the mag- 
netic field is near resonance. 




FIG. 14: Density of atoms in the atomic condensate as pre- 
dicted by our model (solid line) and the Gross-Pitaevskii 
equation (dashed line) during a simulation of a pulse to 
Spuise ~ 156 Gauss. 
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FIG. 15: (Color online) Velocities of molecules during a por- 
tion of a simulation involving a pulse to 156 G. 




FIG. 16: (Color online) Condensed atom density during a 
pulse to -Brea. The inset magnifies the region where the oscil- 
lations in density are most pronounced. 



established in Figure 12 with fewer than 10,000 con- 
densed atoms remaining when iramp = 50 microseconds. 
Figures [THj [TTj and [18] show the number densities of 
condensed atoms, molecules, and noncondensed atoms, 
respectively, for the simulation with Bpuise = -Sres and 
^ramp = 12.5 microscconds. 



as the gradient of the phase. Such instances of very high 
velocities appear periodically, when the molecular field 
declines. We infer that the molecular binding energies are 
transformed to kinetic energy as the molecules dissociate. 

We also performed simulations with Bpuise = -Bros- 
Note that the effective scattering length diverges for this 
value of the magnetic field, and the Gross-Pitaevskii 
equation becomes undefined. The results follow the trend 



V. CONCLUSION 

We have derived a two-channel model of BECs in 
which the constituent atoms interact via a Feshbach reso- 
nance. The model includes first order fluctuations around 
the atomic mean field and is computationally feasible 
in spherical or cylindrical symmetry. We compare our 
model to the results of the collapse experiments at JILA, 
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FIG. 17: (Color online) Molecule density during a pulse to 
Bros- The inset magnifies the same region as the inset in 
Figure [16] 
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FIG. 18: (Color online) Noncondensed atom density during 
a pulse to -Bros- The inset magnifies the same region as the 
inset in Figure [TB) 



which so far always have been difhcult to explain by the- 
ory. 

The model, which accounts for interactions using a 
mean field of bound pairs, approximates the experimen- 
tally measured time to collapse in a simulation of a BEC 
with attractive interactions, but only when the renormal- 
ization incorporates the correct molecular binding energy 
at the cost of incorporating an incorrect mean field en- 
ergy. It is important to note that our model cannot sim- 



ulate the whole regime between an unstably attractive 
gas and a noninteracting gas with a single value of the 
parameter g^. 

The inclusion of three-body effects does not signifi- 
cantly affect the condensate's dynamics leading up to 
collapse, in contrast to recent experimental findings at 
the Australian National University [Tl, which operated 
at much higher densities. Three-body effects may in- 
crease the magnitude of atom loss, especially near and 
after icoUapsc, when the condensate's density is very high. 
Near and after the collapse time, the density is strongly 
localized at the origin, and our simulations cease to be 
numerically tractable due to our choice of a regular grid. 
A finite element approach or other alternative grid may 
allow post-collapse simulations in the future. Such an 
approach could be particularly effective if utilized in the 
center of mass coordinate. 

The model does not quantitatively reproduce the re- 
sults of experiments with a single brief period of repulsive 
interactions. In our simulations of these experiments, a 
significant number of molecules are formed, concurring 
with expectations [3l SI [32] that a nontrivial molecular 
condensate is formed and coexists with the atomic con- 
densate under these circumstances. We found that the 
use of three-body recombination did not improve mat- 
ters. However, the large number of molecules strains 
the assumption in our model that fluctuations around 
the molecular fleld are negligible. Incorporating normal 
and anomalous molecular fluctuations would add consid- 
erably to the computational requirements. Interactions 
of molecules with atoms and with each other, such as 
collisions, may temper the high velocities observed in our 
simulations, possibly resulting in a larger noncondensed 
component and smaller condensed flelds. 
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